Predicting Zonal Flows — A Comprehensive Reynolds-Stress Response-Functional 
from First-Principles-Plasma- Turbulence Computations 
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Turbulence driven zonal flows play an important role in fusion devices since they improve plasma 
confinement by limiting the level of anomalous transport. Current theories mostly focus on flow ex- 
citation but do not self-consistently describe the nearly stationary zonal flow turbulence equilibrium 
state. First-principles two-fluid turbulence studies are used to construct a Reynolds stress response 
functional from observations in turbulent states. This permits, for the first time, a reliable charting 
of zonal flow turbulence equilibria. 



Introduction. — Zonal flows (ZF) in toroidal fusion 
devices are flux-surface averages of layered radial elec- 
tric fields causing poloidal E x B flows with zero poloidal 
and toroidal mode numbers. Stationary ZFs, dominant 
in the core region, are governed by Reynolds stress (RS) 
and reduce the level of anomalous transport by ion- 
temperature-gradient (ITG) turbulence by orders of mag- 
nitude [H-0|. Hence, it is imperative to understand the 
ZF evolution and take their influence into account for 
anomalous transport predictions. Nonlinear analytic ZF 
theories are largely based on wave-kinetics [4-10] and 
wave-kinetic effects have been numerically observed in 
[Tl| . But most numerical studies are restricted to the 
observation of exponential ZF growth and anomalous 
transport reduction 12|, |13| . In order to understand 



the ZF evolution and the radial scale length observed 
in experiments though l4-|l7||. a description for the ZF- 
turbulence equilibrium is necessary, which is not provided 
by contemporary ZF theories. In the following, the time 
evolution of the ZFs is investigated and it is shown that 
the ZF-turbulence equilibrium state is very determinis- 
tic. This indicates that the construction of a determinis- 
tic RS response functional is feasible. The observations 
yield a RS response functional that describes the ZF ex- 
citation, finite saturation and characteristic radial scale 
length and permits, for the first time, a reliable charting 
of ZF-turbulence equilibria. 
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FIG. 1. Top: self-consistent flow pattern (vE y ) = ve v - Mid- 
dle: perp. RS R±. Bottom: shearing rate u = d x {ve v ) • The 
red/black color-coding denotes flows in electron/ion diamag- 
netic drift direction. 



The turbulence equations. — The turbulence is de- 
scribed by the electrostatic two-fluid equations assuming 
adiabatic electrons. The set of equations [3, [l!| for the 
potential, temperature and parallel velocity fluctuations, 
0, Tj and fy with singly charged ions, the minor and ma- 
jor radii r and i?, and equal electron and ion background 
temperatures T e o and T^o is 



D t (0 - (0» - Vx • AVi (20 - (0) + Ti) 

+ d y (/> - e n C (20 - (0) + + e„0||i>ll 
D t (Ti-2(0-(0))/3) 

+ (77* - 2/3) 9 a - 5e„CTi/3 - 2^9^/3 

D t v\\+e v d\\ (20 + Ti) 
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The unit for the coordinates x and y in the radial 
and poloidal directions is the ion gyro radius, pi — 
\/m.iTio/ (eB), whereas the parallel coordinate z = 9 
ranges from — it . , . tt. The parallel length unit is Lu = qR 
(q is safety factor) . The electron adiabaticity relation for 
the density is n = — (0), where the operator (...) de- 
notes a flux surface average, and D t = dt + (z x V^0) • 
Vj_. The unit for the fluctuation quantities e0 and Ti 
is TioPi/L n , the velocity unit is v,ii = c s pi/L n with 
the ion sound speed c s = yj Tio / rm , and the time unit 
to = L n /c s . The density and temperature gradient 
lengths are L n = dr/d (In no) and = dr/d (In T^q). 
The parallel heat conductivity K% is chosen to obtain 
damping rates similar to kinetic phase mixing 
e„ = 2L n /R, rji = L n /L Tll e v = e n / (2g) are dimension- 
less parameters. For circular high aspect ratio geometry, 
the curvature operator is C = cos z d y + sin z d x and the 
parallel derivative is d\\ = d z — sxd v for magnetic shear 
s. The mechanisms of ZF generation and saturation de- 
scribed by this system have been addressed in Ref. [l9| 
including the ZF shearing properties, decreasing radial 
transport, and the appearance of a Dimits shift [21 1. 

Stationary ZFs are defined by the wave-number K y = 
0, the frequency to — and — (0). Integrating the 
flux-surface-average of Eq. ([1]) over x and subtracting 
the flux-surface-average of Eq. (J3j) times —e n /e v cos z and 
Eq. (f2|) yields an equation for the ZF evolution. Therein 
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the parallel velocity time-derivative is replaced using the 
return-flow relation vu = —2qcos zd x {<$>) (required for 
stationary ZFs to cancel the poloidal flow divergence with 
the parallel one), which results in 



d t ( (1 + V cos(z) 2 ) VE v ) = -d x Rt-t 



(4) 



where Rt = R± — 2qRn is the total RS and £ = 
5e„ (sin zc^T^) /3 a finite larmor radius correction which 
is neglected in the following. The perpendicular RS 
R± = (ve v {vdix + weJ) consists of the E x B-velocities 



ve v = d x (j>, ve^ = —dy<j) and the diamagnetic con- 
tribution Vdix = —d y (4>—(4>)+Ti). The parallel RS 
is R\\ = (cos z v\\veS)- Alternatively, if one consid- 
ers ve v — (d x <fi + d x Ti) as the ZF velocity then R± = 
(ve v ve x ) and f = 0. 

Deterministic flows. — A turbulence study with the 
parameters e n = 1.0, q — 1.5, s — 1.0, r\i — 2.4 and a do- 
main of L x x L y x L z — 240pi x 570pi x 2ttL\\ discretized 
over a grid n x x n y x n z = 512 x 1024 x 32 shows a nearly 
stationary poloidal ZF pattern (Fig. [IJ with a character- 
istic scale length that varies only within a small range 
over time. The ZF-turbulence equilibrium scale length is 
thereby different from the scale length during the initial 
ZF excitation phase. States with arbitrary initial flow 
profiles always decay into the characteristic flow pattern 
demonstrating the robustness of the radial scale length in 
the ZF-turbulence equilibrium. Gyro-kinetic turbulence 
studies using the GYRO code Q qualitatively reproduce 
this ZF evolution, which justifies the use of the fluid ap- 
proach in the following. For large domains the flow and 
RS pattern become more deterministic than for smaller 
ones (e.g. L y = 140pi only) because random fluctuations 
are averaged out to a greater extent. This determinism 
permits the following construction of a RS-functional. 

Functional construction from observations. — Com- 
parison of the RS and ZF patterns in Fig. [T] shows 
that Rj_ appears to be proportional to the shearing rate 
u = d x (vEy)- A study, with the parameters e„ = 1.0, 
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FIG. 2. Time traces at x = L x /(4pi). Top: R± (thin) and 
Q (thick). Bottom: R±/Q (solid) and ao_±u (dashed). Com- 
parison shows that R±/Q ~ ao,±u with coefficient ao,±. 

q = 1.4, ,s = 0.8, r\i = 3.1, a domain of L x x L y x L z = 
\22pi x 244^^ x 2-kL\\ and an artificially maintained flow 
( v E y ) ~ sin (2irxpi/ L x ) that oscillates over time (Fig. 



[2]), exhibits variations of at least one order of magnitude 
in the turbulence intensity Q = {ve^Ti) which coincide 
with large deviations of R± from u. Rescaling of R± 
by Q using a constant coefficient ao.± evidently restores 
the proportionality to u . Comparison of the Q and u 
profiles reveals further that Q is not just a function of 
u but follows variations in u with a delay (maxima in Q 
after u — is reached). This implies that Q should be 
regarded as a separate degree of freedom for the RS-ZF 
system. 




FIG. 3. Time-average of R±/Q (solid) over several flow oscil- 
lations compared to the average turbulence intensity gradient 
, y±d x In Q (dashed) with coefficient j±. 

To remove all RS contributions caused by u, R±/Q 
is time-averaged over several complete flow oscillations. 
The residual RS (Fig. |3j) is evidently proportional to the 
time-average of the turbulence intensity gradient d x InQ. 

To examine the wave-number dependence of the stress 
response a full turbulence study with an artificially en- 
forced flow (ve v ) ~ sin (2n x pi(5 + 15x pi/ L x )/ L x ) is 
used (Fig. [4|. The parameters are e n = 1.0, q = 1.5, 
.s = 1.0, rji = 2.4 and a domain of L x x L y x L z — 
lAOpi x 560pi x 2itLh. Apparently the RS decreases with 
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FIG. 4. Radial profile of an artificial flow study with a radially 
changing wave-number. Top: Comparison of R_l/Q (solid) 
with u (dashed). Bottom: Comparison of R±/Q — ao t ±u 
(solid) with a2,±d x u (dashed) with coefficient «2,x- 

the wave-number K x which confirms that the propor- 
tionality to u is indeed wave-length dependent. Since 
the Eqs. ([J)-© are invariant under the transformation 
x,y,z,n,^,Ti,v\i,Rt -> —x,y,—z,—n,—<f>,—Ti,v\\,—Rt 
the structure of additional terms is restricted. The sym- 
metry only allows polynomials of even orders in K x and 
a fit shows R±/Q — «o.±u ~ ct2,±d^.u with a coefficient 

This symmetry and the previous observations hint that 
additional terms for the functional are polynomial in 
K x ,u and c^lnQ as well, thus further simplifying the 
search for them. 
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The ZF amplitude in self-consistent studies is always 
finite which indicates a nonlinear saturation term for the 
functional. To identify the nonlinearity a turbulence 
study with an artificial flow (ve v ) ~ sin (2nxpi/12L x ) 
and a domain of L x x L y x L z = 280pi x 560pi x 2irL\\ 
but otherwise identical parameters to the previous case 
is used (Fig. [5]). The figure reveals that R±/Q satu- 
rates for large shearing rates. The lowest order nonlinear 
term allowed by the symmetry is u 3 and a fit of u 3 to 
R±/Q — cko,J_ u yields the coefficient /3j_. 




FIG. 5. Radial profiles for an artificial flow study with a 
large amplitude. Top: Comparison of R±/Q (solid) with u 
(dashed). Bottom: Comparison of R±/Q — cto,±u (solid) with 
/3j_« 3 (dashed) with coefficient f3±. 

Identical terms also describe the parallel stress behav- 
ior with coefficients {a 2 i, P, 7 }||- Overall, a functional for 
the total RS is obtained, 

Rt =Q (a Q u - fill 3 + a 2 d 2 x u - jd x In Q) , (5) 

with coefficients a 2 i = &2i,±. — 2ga 2 j.||, /3 = /3j_ — 2q/3u, 
7 = 7^ — 2?7||. The coefficients have to be azi,fi > 
to fit the observed linear ZF growth and saturation. In 
addition 7 > was found in all turbulence studies. The 
functional ([5]) reproduces the RS for self-consistent and 
artificial flow patterns rather accurately including the 
equilibrium set up phase. However, actual solutions of 
Q, using the approximation ([5]), a constant Q, which is 
a good approximation in a self-consistent ZF-turbulence 
equilibrium, and periodic boundary conditions, demon- 
strate that any initial state successively evolves towards 
the largest scale length fitting in the system. This clearly 
disagrees with the observation of the characteristic and 
robust ZF scale length (Fig. [T]) . Apparently one ingre- 
dient is missing in the functional ([5]) to self-consistently 
describe the ZF evolution. 

The missing ingredient. — To identify the missing 
contribution, the ZF behavior induced by the functional 
([5]) is analyzed using a mean-field approximation, d x u rs 
3 (it 2 ) x d x u, where (■ ■ -} x denotes a radial average, which 
yields a ZF growth rate 

To = K 2 X (a - 3/3 (u 2 ) x - a 2 K 2 ) (Q) x . (6) 

The region where Tq > is < K x < K x j t with 
K Xt h = — 3/3 (u 2 ) x / y /c7 2 . Outside this region flows 
are damped. K x j t tends towards zero for increasing 



shearing rates (u 2 ) x , explaining the eventual ZF decay 
except for the largest possible scale length. 

In contrast, the self-consistent ZF behavior (Fig. [T]) 
requires that small and large K x be damped while inter- 
mediate K x continue to grow until the system saturates 
at a finite amplitude. To reflect this growth behavior, 
an additional wave-number dependent term is required 
to incorporate the necessary additional root to Tq = 
resulting in 

T = Kl (a - 3/3 (u 2 ) x + a 2 K 2 x - a 4 K±) (Q) x . (7) 

This formula confines the region of growth to a band 
K x j < K x < K x ,h, where < K x j < K Xt h are the roots 
of r = 0, for sufficiently high shearing rates and 014 > 0. 
Since K x> i increases and K x ^ decreases with the shearing 
rate, the system saturates at K x j — ^ — yj a 2 /2a4. 
The corresponding ZF evolution equation is 

d t u = -d 2 x [Q (a u - /3u 3 - a 2 d 2 x u (8) 
-jd x In Q — a±d x u)\ . 

Numerical solutions of Eq. © always yield a stable state 
with a finite amplitude and scale length as required by 
the phenomenology of the RS-ZF-system, corroborating 
the mean field theory. 

Measurement of the K x -term. — Unfortunately, de- 
spite the very deterministic RS pattern, a least-squares 
fit of (J5J) to the turbulence runs still proved to have un- 
acceptably large errors for the coefficient 04. To ver- 
ify the restriction of the ZF-drive to a wave-number 
band required by mean-field theory, a scenario with a 
fixed background shearing-rate u p = const (providing the 
mean-field component (w 2 )^. ~ ( u p) ) an d a small per- 

turbatory shearing-rate u s ~ sin (2irxpiK x ), (u 2 ) — 
1/2 

0.1 (u 2 )^ , with varying K x is studied. In addition an 
optimal filtering technique is employed to measure the 
stress response to u s . 

Writing the time-average of Rt/Q as 



(fr) {x) = E & ^ ^ + E w + n w ( 9 ) 

where the Pi form a set of Ansatz functionals (Pi (x)) = 
(u s , Up, d x In Q) whereas the Nj (x) constitute a set of 
possible error terms (Nj (x)) = (x,x 2 ) and n (x) repre- 
sents random noise with coefficients an,fij. We estimate 
&i with 



*iLt =E (-i) ^ ( c ^) l3 ^ ( 10 ) 

S = ^P 1 (x i )(C- 1 ) y P 1 (x J ) (11) 
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using a minimal variance estimator. The required covari- 
ance matrix is initially estimated by 



v bg + v l,pPl,ihi + v t,NN l}i N ld (12) 
I 



I 



where P^i = Pi (x. t ) /max (Pi), N tii = N t (xi) /max (iVj) 
with u& ff , u;,p, i>;,jv are "a priori" estimates for the vari- 
ances obtained from observations in a large ensemble of 
turbulence studies. The covariance is iteratively refined 
using the variances computed from coefficient estimates 
at different times. 




FIG. 6. ^-dependence of ct\ in an artificial-flow study with 
u — Up + u s for different primary shearing rates u p equal to 
0.1/0.12/0.15 [a.u.] (thick/thin/dashed). Self-consistent ZF 
wave-number is K Xt „. 

Figure |6] shows measurements of ai for several shear- 
ing rates u p . Evidently shearing-rates u s with small K x 
are damped and the region where ot\ > has a lower 
boundary K x j which increases for larger u p until satu- 
ration, a behavior like the one described by T/K x , thus 
validating formula dHJ). 

Discussion — The response functional for the total 
RS is 



Rt = Q (ctQU + (3u 3 + a 2 d x u 
+^d x In Q + aid^u) , 



(13) 



with coefficients a 2 i = a 2 i t ± — 2qa 2i \\, ft = — 2g/3||, 
7 = 7-1- 2<n||, a a > 0, a 2i>0 ,/3,j < 0. The func- 
tional reproduces the features observed in self-consistent 
ZF patterns (Fig. [1]) and describes the excitation, finite 
saturation and robust characteristic radial scale length of 
the ZFs. 

To gain analytical insight on the coefficients it is in- 
structive to use wave-kinetics, in a drift-wave model sys- 
tem (although strictly speaking there is no radial scale 
separation, a prerequisite for wave- kinetic theory [23j j . 
between the ITG-turbulence driving the ZF and the ZFs 
themselves) . This has been carried out for Eq. (JSJ (with- 
out d x InQ) in 24, 25j. The wave-kinetic equation 
is 

d t N k + V fe w • \7 x N k - V^w ■ \7 k N k = C (N k ) , (14) 

where the adiabatic wave-action is N k oc n(k) \<f>\ k with 
a symmetric fc-dependent coefficient n(k). For drift- 



waves (DW) n(k) = ( 
sion u) — loq + kyVE y 



The disper- 



wo [w = k y / (l + k 2 /) for DW] and the Doppler shift 
kyV E y The "collisionar-term C (N k ) = Aw (N k ,o - N k ) 
describes tendency of N k to evolve into a turbulence equi- 
librium state N k .o- R±_ is then given by 



Ri 



(vE*v Ey ) = Jdk 2 rr k (N k ) 



(15) 



where Tt k = k x k y /n(k). 

An expansion of N k with respect the various orders of 
wave interactions u m , K™u = d x u (m, n € N) and d x In Q 



25 1 and terms of even higher order solves Eq. (fl~4|) and 



one can calculate an estimate for R± assuming d t <C Aw: 



RjQ = J2a 2i dl i u+ bu 3 +cd x \nQ 



(16) 



For N k oc Q and d x N k:0 s» {N ky0 /Q)d x Q the coefficients 
ai,CQ,d\ are independent of Q and given by 



a 2i = (d kx [ir k k y v 2 g i]) k /Auj 2i + 1 
b = (kM 7r fe >,/Aw 3 



(17) 
(18) 
(19) 



c = {-K k v gx ) k / Aw 
with (A) k = J dk 2 A(N k fl) /Q. The numerical values for 





Qo,*|ao 


P\b 


Q2,*|a2 


Ct4,* |d4 


y\c 


R± meas. 


2.1 


-3.1 


1.3 


5.1 


-2.1 


R± w.-k. 


2.5 


-3.5 


0.2 


2 ■ 10~ a 


-0.82 


R« meas. 


0.56 


-1.0 


0.84 


2.7 


-0.13 



TABLE I. Measured and wave-kinetically derived coefficients. 

the coefficients for R± are shown in Table HI The signs of 
the wave-kinetic coefficients agree with the measurements 
for R±. However, using only the coefficients ao~4,j_ and 
qt ,i as discussed in [8lllol. 24, 25 1 results in a total stress 



consists of the local frequency 



functional that always leads to a final ZF pattern with 
the largest wave-length admissible by the boundary con- 
ditions (see discussion of Eq. ([5])). The additional term 
represents a limitation of the flow wave-length, which be- 
comes effective at flow amplitudes approaching the self- 
consistent ones but is absent for small amplitudes. It is 
therefore essential to take the higher order contributions 
by R\\ and by d x u to Rt into account, a fact that has 
been largely neglected in all contemporary wave-kinetic 
ZF theories. 

Outlook — The discussed measurement technique 
can, in principle, be applied to obtain the coefficient de- 
pendencies on the plasma parameters or to investigate 
the initial conditions of turbulence triggered transport 
barriers. Further studies of the stress response may re- 
veal additional higher-order wave-numbcr-dependent or 
nonlinear terms that describe metastable ZF states. This 
would allow a reliable charting of ZF-turbulence equi- 
libria opposite to the standard procedure of turbulence 
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parameter scans where the metastability might not be 
identified. 
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